Traveling Waves in 2D Hexagonal Granular Crystal Lattices 
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We describe the dynamic response of a two-dimensional hexagonal packing of uncompressed stainless steel 
spheres excited by localized impulsive loadings. After the initial impact strikes the system, a characteristic 
wave structure emerges and continuously decays as it propagates through the lattice. Using an extension of the 
binary collision approximation (BCA) for one-dimensional chains, we predict its decay rate, which compares 
well with numerical simulations and experimental data. While the hexagonal lattice does not support constant 
speed traveling waves, we provide scaling relations that characterize the power law decay of the wave velocity. 
Lastly, we discuss the effects of weak disorder on the directional amplitude decay rates. 



I. INTRODUCTION 

Granular crystals can be defined as ordered closely packed 
arrays of elastically interacting particles. The fundamen- 
tal wave propagation in one-dimensional (ID) homogeneous 
granular chains has been a topic of numerous investigations 
over the past few years, see e.g. the reviews yJ-Ol and 
references therein. In ID, granular chains support the for- 
mation of highly localized traveling waves, which can be 
described by relatively simple numerical and analytical ap- 
proaches 01 0, ID]. The highly nonlinear response of these 
crystals stems from the Hertzian contact interaction between 
compressed spheres, and the absence of tensile response be- 
tween the grains [6]. In this setting, the presence or absence 
of static precompression has been used to control the intrin- 
sic wave propagation properties of impulsive stresses |Q] @]. 
In addition to their prototypical Fermi-Pasta-Ulam type func- 
tional form and features (such as the existence of supersonic 
traveling waves), granular crystals have been proposed for dif- 
ferent applications including shock and energy absorbing lay- 
ers Ir7l-ll0ll. actuating devices 111 ill , acoustic lenses 01311 . acous- 
tic diodes lfl4ll and sound scramblers 

Although one-dimensional granular chains have been stud- 
ied at considerable length, homogeneous higher dimensional 
crystals have only recently started to be explored in more de- 
tail. Examples of the topics examined include: the dynamic 
load transfer path in disordered two-dimensional systems lfl7l - 
l20ll . the stress- wave anisotropy in centered-square nonlinear 
chains with different materials ll2ll - l23l . and the redistribution 
of energy in a square lattice with one or more interstitials Ir24ll . 
Only a few studies have investigated the dynamic response 
of basic two-dimensional (2D) hexagonal particle arrays 11251 - 
l28tl . focusing primarily on the weakly nonlinear (i.e. com- 
pressed) system response ll29ll30ll . Although the theme of ex- 
citation propagation is of fundamental interest in its own right, 
it may also be of significant practical interest. For instance, 
in applications such as sound absorption, the goal is often to 
find conditions that promote energy dispersion (or dissipation) 
and it is thus of interest to examine such properties of differ- 
ent types of bead arrangements in higher dimensions. In that 
light, we are interested in the dynamic response of impacting 
an initially at rest hexagonal lattice. 

In one spatial dimension, traveling solitary waves are gen- 



erated upon striking one end of an initially at rest chain of 
particles 01 [H]]. The properties of these solitary waves de- 
pend on the composition of the chain, on the geometry and 
material properties of the particles and on the degree of static 
precompression applied to the chain. The two main analyt- 
ical approximations developed in order to study such trav- 
eling waves are based on a quasi-continuum approximation 
01 0], where a partial differential equation is derived from 
the pertinent lattice model, and a binary collision approxi- 
mation (BCA). The BCA relies on the assumption that at a 
given time, a significant portion of the energy of the traveling 
structure is centered over two lattice sites where the resulting 
equations can then be solved exactly 13111 . We should note in 
passing that in addition to these methods and especially in ID, 
the Fourier space, co-traveling frame analysis of ijU enables a 
computation of the exact traveling wave up to a prescribed nu- 
merical tolerance. In 2D configurations, one can study trav- 
eling structures by observing a row of beads along different 
radial directions from the impact bead. For a square pack- 
ing, it has been shown that quasi-one-dimensional traveling 
solitary waves emerge upon striking the lattice HH, and are 
essentially described by the one-dimensional theory. In other 
arrangements, such as hexagonal packings, we argue that such 
quasi-one-dimensional motion is impossible, since each ad- 
jacent row and column will be affected upon being struck, 
regardless of the striking angle. After initially impacting a 
single bead, the energy will gradually be spread over progres- 
sively more and more lattice sites. Since the energy continues 
to spread to an increasing number of beads, the (energy and) 
velocity magnitude will gradually decrease, and thus a perfect 
traveling solitary wave will be impossible to support. This is 
in contrast to the ID (resp. 2D square) situation, where the 
amplitude of the velocity profile remains almost constant, and 
hence, supports a traveling solitary wave. Thus, hexagonally 
packed granular crystals are far better candidates for applica- 
tions such as sound absorption. 

In this paper, we use numerical, theoretical, and experimen- 
tal tools in order to study how energy is spread throughout a 
hexagonally packed granular crystal lattice upon being struck 
at different "strike angles" and being observed at different 
"observation angles". In addition to confirming the absence of 
a traveling solitary wave, we identify scaling power law pat- 
terns for the decay of the velocity as the wave travels through 
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FIG. 1. Schematic diagram of the experimental setup. The solid 
stainless steel sphere locations are represented with grey circles. The 
confining walls are represented by solid black lines. The striker 
sphere impacts the system as shown, with an initial velocity Vstriker 
in the x-direction. The locations of sensor particles are marked with 
a red circle and outlined in black. The number next to each sensor lo- 
cation indicates the different sensor configurations: 1, 2, and 3 (with 
4 sensors present for any given experiment). 



multiple layers of beads within the structure. In a simple spe- 
cial case, we are able to generalize the BCA into a ternary 
collision approximation (TCA), which yields a reasonable nu- 
merical (and approximate semi-analytical) estimate for the ve- 
locity decay. Lastly, additional numerical simulations incor- 
porating weak disorder in the particle lattice are performed 
in order to more realistically simulate and compare with the 
experimental results. 



II. EXPERIMENTAL SETUP 

We performed experiments on a 2D hexagonal array of par- 
ticles consisting of 11 spheres along each edge of the lat- 
tice (see Fig.[TJ). The spheres were stainless steel (type 316), 
with diameter d — 19.05 mm [33]. For subsequent calcula- 
tions and numerical simulations, the particles were assumed 
to have a Young's modulus E = 193 GPa, Poisson's ratio 
v = 0.3, and density p = 8000 %/m 3 iH. The hexagonal 
particle lattice was confined, but not compressed, by walls on 
all 6 sides, with a hole at the impact location along one edge. 
Since we are concerned with only the primary traveling ex- 
citation, a comparatively "soft" material, delrin, was chosen 
to line the confining walls (E — 3.1 GPa, v — 0.35, and 
p = 1400 fc s/m 3 I32IP in order to delay reflections from the 
system edges. To aid in the lattice particle alignment, a slight 
tilt (< 5°) was imposed in the x-direction of the experimental 
setup. 

To excite the system, a striker sphere, identical to the 
spheres composing the 2D hexagonal granular crystal, im- 
pacted the central particle along one edge. The striker sphere 
rolled through a channel down the inclined experimental setup 
and its initial velocity was calculated from the drop height. To 



measure the motion of individual particles within the system 
we used custom fabricated sensor particles that consisted of a 
miniature triaxial accelerometer (PCB 356 AO 1 with sensitiv- 
ity 0.51 mV /^) embedded in a stainless steel sphere. A de- 
tailed description of the sensor particles can be found in ll32ll . 
Based on the data acquisition system used (NI BNC-21 10 and 
NI PCI-6123 with simultaneous sampling rate at 500 ks /s), 
the number of sensors present in the array during a single ex- 
perimental impact was limited to four (with data acquisition in 
both the x- and y-directions). Therefore, to better capture the 
system response, three different sensor configurations were 
used: (1) along the 0° observation direction, (2) perpendicu- 
lar to and near the impact, and (3) perpendicular to and further 
from the impact (see Fig. [TJ. The sensor in the lattice location 
closest to the impact was present in all three sensor configura- 
tions. 

Previous studies on the effects of weak disorder in tightly 
packed 2D granular arrays showed that the exact system re- 
sponse depends on the initial contact lattice, and is generally 
consistent between repeated impacts on an individual con- 
figuration [22}]. In order to capture the variability caused 
by differences in initial contact lattices, 20 different particle 
packings were assembled and experimentally tested (for each 
sensor configuration, 1-3) with an initial striker velocity of 
Striker = 0.40 m /s. Each of the 20 initial assemblies was im- 
pacted repeatedly to calculate average wave front amplitude 
and arrival time values, which could differ due to slight varia- 
tions in impact conditions, such as exact alignment and speed 
of the striker particle, and minor particle rearrangements . Ad- 
ditionally, experiments for 20 different initial particle pack- 
ings with sensor configuration 1 were performed for striker 

velocities, Striker = 0.25 m /s and Striker = 0.70 m /s. 

The impact conditions of the experimental assembly were 
chosen for experimental ease and consistency. In the numer- 
ical simulations and theoretical considerations below, the ini- 
tial velocity Vq was assigned to the n = bead, that is the 
bead impacted by the striker sphere in experiments (see also 
Fig. |2). Numerical simulations, comparing the system re- 
sponse for an array with Striker = Vo and V n= a = Vq, 
showed that the difference was negligible for the studied sys- 
tem. 



III. MODEL 

In order to model the above experimental setup, we con- 
sider an ordered homogeneous lattice of uncompressed spher- 
ical beads arranged in a hexagonal configuration, neglecting, 
at least for the purposes of considerations herein, dissipative 
effects. The Hamiltonian of such a system has the form, 



H = \ M P 2 m,n + V(\dei + q m +2,n ~ 9m,n|) 



+V(\de 2 + q m +i,n+i - q m .n\) 
+V(\de 3 + g TO+ i, n _i - q m .n\) 



(1) 



FIG. 2. (a) Orientation of the hexagonal lattice and possible striking and observation angles. The numeric labels correspond to the counting 
conventions used along the respective observation angles, (b) Labeling convention for the ternary collision approximation 



where q m ,n G is the displacement of the {m, n} bead from 
its initial position, d is the bead diameter, M is the bead mass, 
ei = [2, 0] T , e 2 = [1, %/3] T , e 3 = [1, -V3] T and 

V{r)=Ahd-rfl 2 (2) 
5 

where A is a parameter depending on the elastic properties of 
the material and the geometric characteristics of the beads fll, 



and the bracket is defined by [r]+ = max(0, r). All parame- 
ters can be rescaled to unity using the transformation, 



q m ,n = dq m , n , t = y Adl/2 t (3) 

where the variables with the tilde represent the solutions in 
the physical scaling (i.e., the dimensional ones). We define 

q?n,n \^m.TL-}Vra^ri\ and qm.n Pm,n L^7n,nj ^n,Tn] ■ ^ 

these variables, the equations of motion have the form, 
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where rj — . (r*) 2 + (rf) 2 for j = 1 . . .6 and 



(4) 



n V - 



sin(O) + y m +2,n ~ Vm,n 
sin(7r/3) + y m +i, n +i - y m ,n 



COS(O) + X m+ 2, n ~ Xm,n 

cos(7r/3) + x m+ x, n +i - x m 
cos(-7r/3) - 5G m -i.,„+i + a; m ,„ rf = sin(-7r/3) - y m _i, n +i + y 

COs(O) ^m— 2,n ~t~ X ran 
COs(7r/3) - X m -l,n-l + X m ,n 

cos(-7r/3) + x m+x n _ a - x m 



r\ = sin(O) - y m -2,n + y m ,n 
rf = sin(7r/3) - J/ m -l,n-l + 
r e = sin(-7r/3) + y m+ i, n -\ - Vm, 



We can solve this system numerically using a standard 4th 
order Runge-Kutta method. The delrin-lined walls present in 
experiments were modeled as spheres with an infinite radius. 
The numerical simulation results for the configuration corre- 
sponding to the experimental setup are shown in Fig. [3] Ad- 



ditionally, numerical simulations were performed for a larger 
array consisting of 45-spheres along the bounding hexagonal 
edge, compared to the 11 -sphere per edge system tested in 
experiments. Figure |4] shows the evolution of the wavefront 
shape over the larger 2D hexagonal system. It is interesting 
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FIG. 3. Numerical simulation results showing the wave front shape evolution for the 11 -sphere edge system tested in experiments. The array 
was centrally impacted along the edge by a striker sphere with initial velocity V s triker = 0.40 m /s. The colorbar indicates particle velocity 
magnitude in m /s. 




FIG. 4. Numerical simulation results showing the wave front shape evolution for a larger hexagonal packing (91 particles along the left, long 
edge, and 45 particles along the shorter edges). The array was centrally impacted along the edge by a striker sphere with initial velocity 
V s triker = 0.40 m /s. The black hexagon indicates the wall locations for the size of the experimental setup. The colorbar indicates particle 
velocity magnitude in m /s. 



to observe therein, in addition to the continuous decrease of 
the wavefront speed, its gradual deformation in what appears 
to be a semi-circular shape. We will return to this point later 
in the text (in the discussion after presenting our detailed re- 
sults). Although this model ignores dissipation and disorder 
in the system, it has been found to provide a good comparison 
with the experimental observations, as we also detail below. 



IV. TERNARY COLLISION APPROXIMATION: A 
THEORETICAL APPROACH FOR 30° ANGLE OF 
OBSERVATION 

System (0]l is extremely complex. Indeed, many theoretical 
studies of two-dimensional lattice equations are for simpler, 
single component systems ll34l - [37ll . although traveling waves 
in fully 2D Hamiltonian square lattices of springs have also 
been examined; see e.g. Oal . To help reduce the complex- 
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ity of the system we consider a setting which is amenable to 
a semi-analytical description. More specifically, we turn to 
an analog of the BCA that has been developed for ID chains 
ifJll l39l - l42ll . We first assume that the beads travel approxi- 
mately along fixed angles upon impact, based on symmetry 
considerations. From a theoretical point of view, it is more 
straightforward to develop the theory for a striking and obser- 
vation angle of 30° since the pattern does not alternate along 
this line of observation, in contrast to what is the case for a 
0° striking, see Fig. [2] Considering the time scale where only 
the beads adjacent to the struck bead are affected results in 
a so-called ternary collision approximation (TCA). Actually, 
there are four beads involved in the distribution of the energy 
at each collisional step, yet we use the symmetry to reduce 
the number of degrees of freedom to three. In the renormal- 
ized system, we have the following TCA for the system of the 
right panel of Fig. [2] 

-[x - xifl 2 - 2cos(6»)[cos(6»):ro - x a f_/ 2 



x = 

Xi = [X 
X a = [COS(6»)X - X a f_/ 2 



i3/2 



where xq is the displacement along the impact direction, x\ is 
the displacement along that same direction of the bead adja- 
cent to the impacted bead, and x a is the displacement of the 
bead adjacent to the impacted bead at the angle 9. By sym- 
metry, the bead adjacent along the —9 line will have the same 
contribution as the x a bead. Now define z = xq — x\ and 
y = cos(#)iz:o — x a . We have then 



z = -2\z 



y= 



3/2 r i3/2 

+ -m+ 

12 12 - 3[y] 3 + /2 /2 



(5) 



where we used 9 — 7r/3, which corresponds to a hexagonal 
packing. Although System (O is remarkably simpler than the 
original System ©, it only lends itself to an exact (analyt- 
ically obtainable) solution in the case where y(t) = uz(t), 
where ui £ K. However, note that this prescription requires 
also that the initial conditions satisfy the same scaling. The 
definition of z and y necessitate (from the corresponding ini- 
tial conditions) that ui = 1/2. However, using the above 

ansatz in System (0 yields z = —(2 + w 3 / 2 )[z]+ 2 , as well as 
the "compatibility condition" 2 +w 3/2 = (l/2 + 3/2w 3 / 2 )/w, 
whose sole real solution is oj w 0.382. From the above, it is 
clear that the compatible solution is close to (although not ex- 
actly) satisfying the initial conditions. In that light, we will 
also use for comparison the exact analytical solution of the 
form: 



t = 2 Fi([2/5, 1/2], [7/5], 



2V2 + u; 3 / 2 z 5 / 2 . 
5E 



(6) 



where E in the relevant hypergeometric function expression 
stands for the constant "energy" of the oscillation associated 
with z(t) and can be computed as E = v 2 /2, where v is the 
initial velocity of the bead x . 

As explained in 13 ill , the main purpose of the BCA is to 
offer an approximation of the velocity of the traveling wave 



that propagates through the chain. In order to compute this 
"pulse velocity" one needs to define the time T n the pulse 
resides on bead n. Let t — be the time when the velocities 
of beads n— 1 and n become equal. Then T„ is the time when 
the velocities of beads n and n + 1 become equal and the pulse 
velocity observed at the nth bead is c n = 1/T n . 

In order to approximate the residence time T n , the BCA 
combines two points of view. First, it is assumed that the 
interaction is not instantaneous, such that the velocity of the 
impacted bead will gradually decrease. However, the BCA 
is only valid if two beads are involved in the dynamics, 
which is clearly not the case over the entire residence time 
T n . Thus, strictly in terms of the BCA, it is not clear at all 
how to initialize the next step in the iteration. To circum- 
vent this issue, it is supposed that the interaction is instan- 
taneous such that the conservation of kinetic energy and mo- 
mentum can be used to compute the velocity that bead n will 
emerge with after the collision with bead n — 1. For exam- 
ple, for a heterogeneous ID chain, these conservations yield 
V n = 2T4.-i/(l + mn/mn-i), where m n is the mass of bead 
n lHH. Applying this relationship recursively yields 



Vn 



+ rrij/mj—i 



(7) 



Thus, the reduced equations in the BCA setting at the nth step 
are initialized with z(0) = 0, z(0) = V n . Although it appears 
somewhat awkward to adopt these two viewpoints in order to 
make the approximation work, the results in ID chains have 
been surprisingly good Il3ll 139144211 . For this reason we carry 
out a similar procedure to compute the pulse velocity of the 
traveling wave in the hexagonal lattice. 

Due to the complexity of the hexagonal system there is no 
clear way to obtain an analog of relationship (|7|i using con- 
servation of kinetic energy and momentum. However, in both 
the ID BCA and 2D TCA setting, we can compute the time it 
takes for the initially impacted bead to obtain a velocity equal 
to its adjacent partner. For a homogeneous ID chain, the BCA 
predicts this to occur for a velocity equal to 50% of the initial 
velocity. Thus, to understand what contribution the additional 
beads in the hexagonal configuration will absorb, we use the 
reduced equations (|5]l to compute when the velocities of beads 
and 1 are equal (see Fig.|2]for the labeling convention). The 
ratio of the velocity at this time compared to the ID case is one 
way to quantify the fraction of momentum that is absorbed by 
the additional beads. Call this value c. For example, suppose 
the velocity at the time of intersection was predicted to be 0.45 
in the TCA and 0.5 in the BCA. Then c = 0.9, since the bead 
only reached 90% of the value it would have in the absence of 
the additional beads. If in addition, we assume that this ratio 
remains constant as the traveling structure propagates through 
the lattice, we expect that bead n emerges with the velocity 
V n = cV n -\ after collision with bead n — 1. Applying this 
relationship recursively yields, 

V n = V c n , (8) 

where Vq is the impact velocity. Equation ([8]) will now play 
the role of (0 for the hexagonal system. In the rescaled system 
we found c w 0.8452. 
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Given the repetition of the fundamental building block of 
the right panel of Fig. [2] we now simply solve the TCA Sys- 
tem (j5]l with initial values defined through (HJ and compute 
the intersection time. Doing this for several bead locations 
yields a discrete set of approximations for the pulse velocity 
observed at those bead locations. See the left panel of Fig. [5] 
for an example with Vq = 0.1. There, it can be seen that 
such an approximation of this average pulse velocity favor- 
ably compares with the full numerical results for the hexago- 
nal lattice. 

To amend the TCA to a 0° strike (but still observing along 
the 30° line) we need to understand how the energy is trans- 
ferred among the first two beads that are in contact with the 
impact bead (see Fig. |2). Due to symmetry considerations, it 
is reasonable to conjecture that the velocity contribution along 
the 30° (resp. -30°) will be half of the impact velocity. This 
would be the case if kinetic energy and momentum were con- 
served. This relation was observed to be fairly accurate (see 
the right panel of Fig. |7]l. At the 0° strike and 30° observation 
angle we have experimental data available for comparison, see 
the right panel of Fig. [5] Due to the limited number of sensors 
available, we simply present the arrival time data (which is the 
sum of the relevant residence times). We should also note that 
the experiments use the arrival time at the first sensor position 
in Fig. [TJto define t = 0. Thus, in order to compare the ex- 
perimental values to the TCA and numerical simulations we 
estimate the arrival time at the first sensor using the numerical 
simulation. 



V. A NUMERICAL & EXPERIMENTAL STUDY: 
VELOCITY DECAY RATE FOR VARIABLE 
STRIKING/OBSERVATION ANGLE COMBINATIONS 

In applications, not only is the time of arrival of the trav- 
eling wave important, but the magnitude of its amplitude as 
well. The TCA, although simple enough to afford analytical 
approximations, ignores the (weak) dependence of the decay 
rate c on Vq and is restricted to an observation angle of 30°. 
We wish to have a more accurate and robust description of 
the velocity decay rate and thus we now turn to a (strictly) 
numerical study of other combinations of striking/observation 
angles. We corroborate the obtained qualitative and quantita- 
tive results by means of experimental observations. 



A. 30° Strike: 

In panel (a) of Fig.|6j the magnitude of the velocity of each 
bead along the 30° line is shown against time. Notice how the 
maximum velocity obtained for each bead decreases as we 
move further from the impact point, as expected. We fit these 
maximum points with a function of the form Vic™ -1 where c 
is a decay parameter that will now depend on the initial strike 
velocity Vq (in contrast to the situation above which had c in- 
dependent of Vq). We start the fit at V% since the initial impact 
will only have a velocity component, and thus the dynamics of 
the collision of the n — and n = 1 bead are inherently dif- 



ferent than the dynamics of the n and n+1 beads for n > 0. In 
the latter case, the bead will have a velocity and position com- 
ponent. In panel (b) of Fig.|6]we fit the decay rates c = c(Vq) 
with a linear function mVo + b with two free parameters m, b. 
We found that m ps 0.257 and b ps 0.832. In panel (c) of 
Fig. [6] we obtain the relationship between Vq and V%, which 
has the form V x = aV 2 + /3V Q . We found a w 0.206 and 
f3 ps 0.629, such that the decay between Vq and V\ is greater 
than that between V n and V n +i for n > 0. We have then that, 

V n = V {aV +p){mV +b) n -\ n>0 (9) 

where V n is maximum velocity of the bead lying n beads away 
from the strike point on the 30° line. Clearly this description 
of the velocity decay is more accurate than © since the de- 
pendence on the initial impact velocity is taken into account. 
Nevertheless, we still observe that a power law decay of the 
maximum velocity of each bead provides an accurate descrip- 
tion of the dynamics along the 30° observation angle in the 
hexagonal chain. 

When observing along the zero degree line it turns out 
that the fitting is more accurate when starting at bead 2, see 
Fig. [2] for the labeling convention. Panel (a) of Fig. Q shows 
a relavent example where the velocity fitting starts at bead 2, 
but for a 0° strike (which is analyzed in Sec. IV Bt . Since the 
relationship between V\ and Vq was linear, it is natural to sup- 
pose that the relationship between V2 and Vq is quadratic, i.e. 
V2 ~ olVq + (3Vq + 7 (see panel (d) of Fig. |6]for example). 
We have then that, 

V n = Vq (aV 2 + /3V + 7 )(mVo + &)"~ 2 , n > 1 (10) 

where V n is the magnitude of the maximum velocity of the 
nth bead from the strike point on the 0° line. Our fitting 
yielded m ps .3039, b ps .8756 and a ps -0.4788, (3 ps 
-0.0256, 7 ps 0.1936. 

B. 0° Strike: 

Finally, we carry out the procedure for the case of striking 
at a 0° angle. As mentioned above for a zero degree obser- 
vation, the first bead behaves somewhat differently (see panel 
(a) of Fig. |7), and therefore we start the fitting again at bead 
2. We found that our scaling relations continue to be valid 
for a 0° strike, but now with m ps 0.1766, b ps 0.8878 and 
a ps -0.2643, ,5 pa -.0969,7 ~ -2152. The optimal scaling 
along the 0° observation angle has the same form as (ITdb . 

When observing at a 30° angle, the formula that yielded the 
best fit had the form, 

V n = V {aV 2 + PVq + 7 )(mVo + &) n ~\ n > (11) 

with m ps 0.1484, b ps 0.8393 and a ss -0.1182,/? ps 
0.0203, 7 « 0.5192. However, for the range of initial veloci- 
ties used, the relationship between V\ and Vq is approximately 
V\ = Vq/2 (see panel (b) of Fig.[7j which yields the simplified 
formula, 

V n = y(mVb + &) n -\ n>0. (12) 
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FIG. 5. (a) Observed pulse velocity of the propagating structure in the ideal parameter case at various bead numbers n for a 30° strike and 30° 
observation angle. The pulse velocity is defined as c n — 1 /T n where T n is the nth residence time. The blue dots represent the calculated time 
based on numerical simulations and the red circles are the prediction based on the TCA. The black triangles denote the exact analytical solution 
to the TCA system (which, however, is only approximate). The initial strike velocity was Vo = 0.1. (b) The arrival time of the pulse for the 
experimental parameter values at a 0° strike and 30° observation angle with Vo = 0.4 m /s. The blue dots represent numerical simulations, the 
red circles are the TCA, the black triangles denote the exact analytical solution to the TCA, and the points with error bars at n — 6 and n = 10 
are the experimental values. 




FIG. 6. Illustration of the procedure to obtain the scaling law for the decay of the velocity profile. This example is for strike and observation 
on the 30° line, (a) Magnitude of velocity of beads — 12 versus time in the ideal parameter case with Vo = 0.05. The maxima are fit with 
a function of the form Vic™ -1 , where Vi is the maximum velocity of the first bead (labeled in the plot). This processes is repeated for several 
values of Vo to obtain a sequence of decay rates c(Vb). (b) The sequence of decay rates c(Vb) are fairly accurately fit with a linear function 
(see the text), (c) The relationship between Vo and Vi is computed, and can be either linear (like in this case) or quadratic. Finally, putting 
these relationships together yields a formula for the maximum velocity at a given bead, see e.g. Eq. l|9). (d) For a 0° observation, the fitting 
begins at bead 2, where the relationship between the initial striking velocity Vo and Vi is quadratic. 
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FIG. 7. (a) magnitude of the velocity of beads — 12 along the 0° line upon being struck at a 0° angle. In this case, the fitting is more 
accurate starting at the second peak V2. (b) The relationship between the initial striking velocity Vo and Vi for a 0° strike and 30° observation 
is quadratic. For this particular case, we see that Vi ~ Vo /2. 



Note, in order to obtain the scaling law for arbitrary param- 
eter values, one simply replaces each occurrence of V with 

Vy/MjA/d 5 / 4 . 



VI. NUMERICAL & EXPERIMENTAL STUDY 
INCORPORATING WEAK DISORDER 

A. Preparation Scheme 

Numerical simulations incorporating weak disorder were 
also performed for the 2D hexagonal system described in 
Section HIl Similar to ll22ll . weak disorder was incorporated 
in the simulations by assigning each particle in the array a 
random diameter, based on a normal distribution with mean 
fi — d = 19.05 mm and standard deviation a = 6tol, where 
tol = 0.0127 mm is the manufacturer specified diameter tol- 
erance of the particles used in experiments. The initial rest- 
ing positions of the particles in the weakly disordered lattice 
were found in two steps. First, the particles were given an 
initial hexagonal lattice spacing assuming d = d + tol, to 
avoid large repulsive forces between overlapping particles. To 
bring the particles in contact, a 5N force was applied to each 
of the edge particles and an artificial damping (of the type 
introduced in i43lo was used to settle the random particle mo- 
tion. Secondly, the wall positions were found based on the 
slightly compressed settled configuration. Then, gravity was 
introduced along the x-direction (in agreement with the ex- 
perimental tilt), the 5N force was removed, and the particles 
were again settled through a similar, but weaker damping pro- 
cess. Once the initial positions were obtained, the settled array 
was impacted with a striker sphere, Striker = 0.25, 0.40, or 
0.70 m /s, along one edge of the system. This process was 
performed for 20 different initial lattice configurations. For 
the subsequent discussion, we will use the terms "ideal" and 
"weakly disordered" to distinguish between numerical simu- 
lations where the spheres are all assigned the same diameter, 
and those just described with a variable diameter. 



B. Discussion of Results 

Based on the numerical simulations of the ideal hexago- 
nal lattice (Figs. [3] and |4), the wavefront shape can be ini- 
tially described by a hexagonal pattern which gradually tran- 
sitions into a more circular shape. This transition occurs as 
the scale of the structure changes from being comparable to 
the lattice "spacing" to being much wider than that. It is 
worthwhile to mention that this is an intriguing feature that 
perhaps a quasi-continuum or homogenization type approach 
(see e.g. lfl2TI and references therein for recent work on the 
subject) may capture and is left as an interesting open prob- 
lem. Figure |8] compares the initial wavefront shape observed 
in the ideal numerical simulations with those observed in the 
experiments and numerical simulations incorporating weak 
disorder. The average wave front arrival times from both ex- 
periments and numerical simulations incorporating weak dis- 
order clearly show the initial hexagonal wave front propagat- 
ing though the structure. The average arrival times in ex- 
periments are slightly longer than those predicted from the 
ideal numerical simulations. Physically, this can be under- 
stood since the dissipation present in experiments results in a 
decreased wave amplitude, consequently decreasing the wave 
front speed. Additionally, imperfections in the experimental 
contact lattices also result in redirection of wave amplitude 
and delays in the signal arrival time 1E2I1 . As Fig. ^indi- 
cates, the relative arrival times in the weakly disordered simu- 
lations are notably longer than both the predicted values from 
the ideal simulations as well as the experiments. This is a bit 
surprising considering that the simulations do not incorporate 
dissipation, thus the time delay is entirely due to the effects 
of the particle misalignments. The discrepancy between the 
arrival times in experiments and weakly disordered numerical 
simulations suggests that the tolerance values used to simulate 
imperfections in the contact lattice are fairly conservative. 

Figure [9] compares the wave front amplitude decay along 
the 0° and —30° observation directions after a 0° strike angle 
for experiments, ideal numerical simulations, weakly disor- 
dered numerical simulations, and the corresponding Eqs. (|T0l 
and (fTTT i. Overall the trends observed for all approaches are 
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FIG. 8. Wave front shape based on the relative wave front arrival time with respect to the particle marked with a white x. (a) Numerical 
simulation of ideal system (without disorder or dissipation), (b) Average wave front arrival time calculated from experiments, (c) Average 
wave front arrival time calculated from simulations with weak disorder. Black circles in the left and right panels indicate experimental sensor 
locations. The colorbar indicates the relative wave front arrival time in milliseconds. 



in good agreement. The presence of weak disorder in the nu- 
merical simulations results in decreased amplitude transmis- 
sion along the 0° observation direction and increased ampli- 
tude transmission along the —30° observation direction, com- 
pared to the ideal simulations. This altered distribution of 
wave front amplitude can be seen in Fig. [9] where the mean 
amplitude from simulations incorporating weak disorder falls 
above the ideal simulation values along the ±30° observation 
direction and below the ideal simulation values along the 0° 
observation direction. The 30° observation direction repre- 
sents a line of spheres directly in contact, while the wave front 
must travel though a zig-zag of particle contacts along the 0° 
direction. Therefore, the number of contacts, or possible am- 
plitude scattering points, is greater along the 0° observation 
direction, which could help physically explain this amplitude 
redistribution phenomenon. Although the presence of dissipa- 
tion in experiments results in lower average amplitudes com- 
pared to the weakly disordered numerical simulations, which 
neglect dissipation, we also observe this trend in average ex- 
perimental amplitude values (see the middle panel of Fig. [9] 
for V = 0.40 m/ a ). 

VII. CONCLUSIONS AND FUTURE CHALLENGES 

We presented a systematic study of the dynamic response 
of a 2D hexagonal, highly nonlinear lattice excited by an 
impulse. In this 2D hexagonal setting, because of the ever- 



increasing number of neighbors over which the energy is dis- 
tributed, no genuine traveling wave excitations, i.e. with con- 
stant velocity, have been found to persist. The propagating 
pulse energy has been found to decay as a power law, both in 
our numerical computations and in our experimental observa- 
tions. Detailed expressions were provided to describe these 
power laws at different angles of strike and of observation. In 
a special case (of 30° strike and 30° observation, according 
to our presented classification), a generalization of the binary 
collision approximation (dubbed the ternary collision approx- 
imation) was presented and utilized to give simple numerical 
and even approximate analytical expressions for the bead evo- 
lution. Lastly, the effects of weak disorder on the propagating 
wave structure were examined; the average spatial amplitude 
values from numerical simulations incorporating weak disor- 
der were in good agreement with experimental values as well 
as the corresponding fitted equations derived from numeri- 
cal simulations on the ideal hexagonal granular crystal. This 
agreement reveals that the level of disorder present in experi- 
ments does not cause significant deviation of the propagating 
wave structure from the predicted system response. 

While this is a fundamental step towards an improved un- 
derstanding of non-square lattices, there are numerous addi- 
tional tasks to be considered. While invariant traveling so- 
lutions may not exist, the possibility of existence of self- 
similar decaying solutions of the discrete model or perhaps of 
a quasi-continuum approximation thereof cannot be excluded 
and should be further considered. Additionally, while here we 
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FIG. 9. Maximum particle velocity with respect to sphere position denoted by n. Panels (a), (c), (e): Strike and Observation. Panels (b), 
(d), (f): 0° Strike and -30° Observation. Initial velocity: (a), (b) Vo = 0.25 (c), (d) Vo = 0.40 and (e), (f) Vo = 0.70 Shaded 
grey regions indicate ±1 standard deviation, a, and dotted grey lines indicate the average value, n, from the 20 simulations incorporating weak 
disorder. The solid red lines represent Equations [TJj]and [TTJ The black dots and error bars represent the mean and standard deviation from 
experiments (Vo = 0.25, 0.40, and 0.70 m / s for 0° Strike 0° Observation and Vo = 0.40 ™/s for 0° Strike 30° Observation). 
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have concerned ourselves with the highly nonlinear limit of 
the model (i.e., no precompression), there has recently been 
a surge of activity in connection to linear or weakly nonlin- 
ear properties of non-square (such as honeycomb or hexag- 
onal) lattices. This is due to some of the remarkable linear 
properties of these lattices, such as, e.g., conical diffraction 
and the existence of Dirac (diabolical) points examined in 
both the physical B44 14511 and mathematical [46] communi- 
ties. It would be particularly relevant to consider the possi- 
bility of such spectral and nonlinear features in the precom- 
pressed variant of the present chain. Some of these possibil- 
ities are currently under investigation and will be reported in 
future publications. 
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